AD-A079  26H  MARYLAND  UNIV  C0LLE6E  PARK  INST  FOR  PHYSICAL  SCIENCE— ETC  F/6  12/1 

RELIABLE  ERROR  ESTIMATION  AND  MESH  ADAPTATION  FOR  THE  FINITE  EL-ETC(U) 
APR  79  I  BABUSKA  »  M  C  RHEINBOLDT  N000U-77-C-0623 

UNCLASSIFIED  BBN-910  nl 


I  »»  I 

^079264 


Institute 


"SdiNoloqy 


Technical  Note  BN-910 


Nrmrsity  of 

I  I  larylond 
College  Fbrk 


tation 


ite  Element  Method 


I .  Babuska 


W.  C.  Rheinboldt 


APPROVED  FOP  PUBLIC  KFL2JLS* 
DIS1RIBUIIQII  UZiLIkllTliD 


April  1979 


4'  . 

Rel  1 


able  Error  Estimation 


and 


Mesh  Adaptation 


fr 


for  the  Finite  Element  Method  , 


by 


Werner  C./ Rheinboldt 


:  1)  T-e.  tk  hi  tojl  /vot*, 


\l )  IjtL 


1)  This  work  was  in  part  supported  under  DOE. Contract  E(401)3443',  NSF  Grant 
MCS-78-05299,  and  ONR  Contra^N0$l  4-77-0-^623^  £ 

Institute  for  Physical  ScierTce^nd  Technology^Uni versity  of  Maryland,  '** 
College  Park,  MD  20742 


2) 


3)  Department  of  Mathematics  and  Statistics,  University  of  Pittsburgh, 
Pittsburgh,  PA  15260 


1.  Introduction 


In  the  past  two  decades  computational  structural  analysis  and  design 
have  developed  into  very  extensive  subject  areas.  A  major  contributing 
factor  certainly  is  the  evolution  of  the  finite  element  method  into  a 


powerful  tool  for  the  solution  of  a  broad  range  of  problems..  In  assessing 
trends  for  the  future  of  these  areas  (see  eg  [1],  [2],  [3])  w£  need  to 
take  account  of  some  of  the  questions  for  which  there  are  as  yet  no 


satisfactory  answers. 

w^In  this  article  we  address  two  such  questions  which  may  be  broadly 


phrased  as  follows:  ^ 

(a)  How  can  we  calculate,  at  reasonable  cost,  reliable  estimates 
of  the  error  of  computed  finite  element  solutions?  », 

^  ^  LMb)  How  should  we  set  up  the  overall  solution  process  so  as 
to  achieve,  with  reasonable  certainty  and  efficiency, 
a  solution  for  which  the  accuracy  either  falls  into  a 
prescribed  range  or  is  best  possible  for  the  allotted 
computational  cost?  , — 

r- 

In  this  form  both  questions  need  further  elaboration. 


In  practice,  the  accuracy  of  a  finite  element  solution  is  usually 
measured  in  relation  to  the  exact  solution  of  the  underlying  mathematical 
problem,  as  for  instance,  an  elasticity  problem.  In  many  instances  there 
is  interest  also  in  the  error  of  the  mathematical  model  itself.  In  [1]  we 
presented  seme  comments  on  the  latter  question,  but  here  we  shall  concentrate 
only  on  the  first  type  of  error. 

Any  measure  of  the  accuracy  involves  the  choice  of  a  norm.  Different 
norms  show  up  different  aspects  of  the  error  and  in  practice  it  is  important 
to  allow  for  a  choice  from  among  any  reasonable  collection  of  norms.  This 
represents  a  demanding  requirement  on  the  design  of  the  error  estimates. 


In  practical  engineering  calculations  the  required  accuracy  usually 
is  not  very  high;  in  fact,  a  range  of  5-10%  is  often  typical.  In  that 
case,  the  efficiency  index  0=  error  bound/exact  error  should  satisfy,  say, 

1.0  <0  <  1.3  for  the  bound  to  be  realistic.  For  accuracies  beyond,  10% 
the  value  of  0  may  become  worse  but,  on  the  other  hand,  0  should  tend  to 
one  when  the  accuracy  becomes  better  and  better.  This  points  to  the  possibility 
of  relaxing  the  requirement  for  a  guaranteed  bound.  In  fact,  we  may  consider 
accepting  estimates  e  of  the  exact  error  ||e||  for  which 

e/||e||=l+0(||e||a)  as  |  |e|  |  -*■  0.  Here  the  bound  in  the  0-term  need 
not  be  explicitly  known  but  the  exponent  a  >  0  should  be  reasonable;  for 
instance,  a  =  2  would  be  rather  desirable. 

A-priori  error  bounds,  as  important  as  they  are  for  the  theory,  are 
rarely  very  realistic  for  practical  problems.  Thus  we  have  to  consider 
a-posteriori  estimates  which  utilize  information  obtained  during  the  solution 
process.  At  the  same  time,  such  estimates  open  up  the  possibility  of  an 
adaptive  control  of  the  solution  process  itself  to  achieve  the  goal  contained 
in  the  second  question  (1.1). 

Generally,  any  overall  assessment  of  the  efficiency  of  the  finite 
element  solution  of  a  practical  problem  should  include  not  only  the  direct 
computational  costs  but  also  the  expenses  for  data  preparation,  etc.  A 
particularly  costly  part  is  here  the  construction  of  the  finite  element  meshes 
which  should  be  automatized  as  much  as  possible.  Moreover,  a  solution  with 
optimal  accuracy  within  a  given  range  of  computational  effort  can  hardly 
be  achieved  without  the  use  of  a  responsive  adaptive  mesh-design  algorithm. 

Any  efficient  answer  to  the  second  question  (1.1)  can  only  be  achieved 
by  a  complete  integration  of  the  error  estimates  and  adaptive  controls  into 


the  software.  For  this,  special  attention  has  to  be  paid  to  data  management 
and  the  overall  program  design.  Here  modern  developments  in  computer 
science  have  to  be  taken  into  account.  Thus  a  satisfactory  answer  to  (1.1 )(b) 
will  have  to  integrate  results  from  engineering,  mathematics,  and  computer 
science. 

In  this  paper  we  present  some  answers  to  the  two  questions  (1.1).  In 
[4]-[6]  a  mathematical  theory  of  a-posteriori  error  estimates  for  finite 
element  solutions  has  been  given.  Instead  of  discussing  this  theory  in 
general  we  present  here  its  principal  aspects  in  the  case  of  several  model 
problems.  More  specifically,  in  Chapters  2  and  3  we  consider  certain 
one-dimensional  and  two-dimensional  linear  problems,  respectively,  and  present 
a  number  of  experimental  results  which  show  the  effectivity  of  the  bounds. 

From  a  theoretical  viewpoint  there  are  still  various  open  questions  which, 
in  Chapter  3, are  formulated  as  a  list  of  conjectures.  In  Chapter  4  we 


summarize  our  current  work  on  an  experimental  finite  element  solver  which 
is  aimed  at  answering  question  (1.1 )(b).  Finally,  in  Chapter  5  we  give  an 
extension  of  the  results  to  nonlinear  problems  and  of  their  integration  with 
modern  continuation  techniques.  Experimental  results  for  a  particular 
nonlinear  two-point  boundary  value  problem  show  that  the  combination  of 
continuation,  error  estimation,  and  adaptive  mesh  construction  should  be  an 
extremely  effective  approach  to  the  computational  solution  of  many  practical 
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2.  One-Dimensional  Linear  Problems 

2.1  The  Model  Problem 

On  the  unit  interval  I  =  [0,1]  c  we  consider  the  elliptic  boundary 
value  problem 


(2.1) 


d_ 

dx 


b(x)u  =  f ( x )  ,  x  e  I  , 


u(0)  =  u(l )  =  0 


Specific  smoothness  conditions  for  a,b,f  will  be  given  later;  but  we  require 
always  that 

(2.2)  a  =  sup  <  «  ,  a  =  inf  a(x)  >0  ,  *  >  b  >_  b(x)  ^  0,  Vx  e  I 

xel  xel 

As  indicated  in  the  Introduction,  the  problem  formulation  is  incomplete 
without  a  specification  of  the  norm  to  be  used  in  measuring  the  desired 
accuracy.  As  an  example,  we  consider  the  Lp-stress -energy  norms- 

(2.3)  | |u| | E,p  =  |  |  [(a(u')2)1/clpdx  |  ,  for  1  s  p  <  « 

(2.4)  |  |u||r  „  =  sup  i  (a(u‘ )2)1/2  . 

L’  xel  1 

Let  uQ  denote  the  exact  solution  of  (2.1),  (2.2)  and  u(a)  any  approximate 
solution  obtained  by  the  finite  element  method  for  a  given  mesh  a  on  I. 


Then  we  have  the  following  tasks: 


(a)  Compute  reliable,  close  bounds  of  the  error  ||un  -  u(a)[L 
under  the  norm  (2.3),  (2.4)  with  a  specified  value  of  p. 

(b)  Design  an  effective,  adaptive  algorithm  to  compute  u(a)  so  that 
ei ther 

(2-5)  ||Uo  -  uU)]|£-p  <ri!u0liE>p 

for  a  given  tolerance  x  >  0,  or  that  ||un  -  u(a)|L  is 

U  t  «  P 

minimal  for  a  given  total  computational  cost. 

2.2  Error  Indicators  and  Estimators  -  I 

Let  P(<)  denote  a  family  of  partitions  A  of  I, 

A  A  A 

A:  0  =  *o  <  *1  <  ' ‘ '  <  xm  =  m  =  mU). 

(2.6)  Ij  =  [xj_1 ,  Xj]  ,  hj  =  Xj  -  Xj_1  ,  j=l,...,  m( a) 

FT(a)  =  max  h^  ,  h(A)  =  min  hA  , 

j  J  j  J 

which  satisfy  the  ^-regularity  condition 

(2.7)  h(A)  _>  c  F(A)K  ,  Va  e  P(x), 

o, 

with  a  constant  c  >  0  independent  of  A.  For  each  Ac  P(<)  let  S  be 
the  set  of  all  continuous  functions  on  I  which  are  linear  on  each  subinterval 
Ij  ,  j=l,...,  m(a),  and  zero  at  the  endpoints  x^  =  0  ,  x^  =  1.  Then, 
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as  usual,  the  finite  element  element  solution  u(a)  e  is  defined  by  the 
condition 


(2.8) 


[a  u(a)'  v'  +  b  u(a) v  ]dx  = 


<•  1 


fv  dx 


Vv 


We  use  a  local  analysis  to  construct  estimates  of  the  error  e(i)  *  Uq-u(4) 

under  any  specified  norm  ||.|j_  •  More  specifically,  from  locally  available 

^  »P 

A 

information  we  compute  an  error  indicator  n-  for  each  mesh-interval 

*3 

Ij  »  j=l,...,  m(a).  These  indicators  are  then  combined  appropriately  to 
yield  an  error  estimator  ^(a)  which,  in  some  sense,  represents  an  approx¬ 
imation  of  the  error: 


(2.9) 


!e(A)|lE>p  -  e(A). 


Finally,  again  on  the  basis  of  the  locally  available  data,  we  calculate  an 
error  corrector  ^(a)  such  that 


(2.10) 


e ( A )  |  |E^p  <  e(A)  +  4>(A) 


Note  that  in  the  evaluation  of  the  quantities  (n->  ,  e(A)  ,  *(a)  only 

J 

information  is  to  be  used  that  is  readily  computable  during  the  computation 
of  u(a)  itself. 

The  estimators  and  correctors  are  said  to  be  asymptotically  exact  if, 
under  reasonable  assumptions  about  the  solution  and  the  partitions,  we  have 


e  (a 
"ilAl 


E.P 


1  , 


TTeT 


E.P 


(2.11) 


0  ,  as  h(  a)  -*•  0,  A  c  P(< ) . 


Furthermore,  we  call  e(a)  an  upper-estimator  ,  or  lower-estimator  if 

(2.12)  | I e( A) | L  <  K  e( a)  ,  Va  e  P(<), 

or 

1 1 e( a) I  I p  >  K  e(a)  ,  Va  e  P(<), 

i 

respectively,  with  constants  K,K  which  are  independent  of  A  and  u(a),  but 
which  may  depend  on  other  input  data.  Thus  7(a)  =  e(a)  +  <j>( a)  always 
represents  an  upper-estimator. 

The  quotient 


(2.14) 


is  called  the  efficiency  index  of  the  estimator  e(a).  As  outlined  in  the 
Introduction,  we  want  0  to  be  as  near  to  one  as  possible.  Since,  in 
practice,  0  is  not  available  this  means  that  we  should  analyse  general 
conditions  under  which  0  approximates  one. 


2.3  Error  Indicators  and  Estimators  -  II 

For  the  problem  (2.1),  (2.2)  we  assume  now  that  a',  b,  and  f  are  p-th- 
power  integrable.  Then  the  residuals 

(2.15)  r.(x)  =  f(x)  +  a' (x)u(A)' (x)  -  b(x)u(A)(x),  Vx  e  1^  ,  j=l,2,...,  m, 

J  J 


(2.16) 


aj  ■  *<R,  + 

* 

r.  =  r^(x)(x-x^_1 )(Xj  -x)dx 


we  introduce  the  following  two  types  of  error  indicators  for  the  norm  (2.3) 
with  some  specified  p,  1  <  p  <  <*>  ; 


(2.17)  ndU)-£<P*l)',/p  a:'/2  hj  [  |  |rjtx)  |pdxl1/p  . 

’j 

>  j  =  l  ,  m( a) 

(2.18)  n . (a)  =  3(p+l)'1/p  aT1/2  hT2  +  1/p  \7.\ 

J  J  J  J  / 


The  corresponding  error  estimators  are  defined  by 

m(  a) 

(2.19)  e( A)  =  [  £  ni(A)P11/p  ,  1  <  P  <  - 

J«1 
m(  a) 

(2.20)  c(  A)  =  [  Y.  n .  (a)P1  /P  ,  1<P<- 

J-l  J 

It  is  readily  seen  how  to  define  analogous  quantities  for  the  case  p=®. 
Then  the  following  result  holds: 


Theorem  2.1:  If  a',b,f  t  Lp(I)  then  e( a)  is  an  upper-estimator  and  c(a)  a 
lower-estimator. 
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Suppose  that  the  smoothness  conditions  for  a,b,f  are  strengthened  to 

(i)  a  c  C2(I),  b,f  e  C1 ( I ) 

(2.21) 

(ii)  Uq'  has  only  finitely  many  simple  roots  in  I; 
then  theorem  1  can  be  improved  as  follows: 

Theorem  2.2:  Under  the  smoothness  conditions  (2.21)  and  for  a  e  P(  <)  the 
estimators  e(a)  and  e(a)  are  asymptotically  exact. 

With 


(2.22) 


^  3  inf  sup  ja(x)  -  t|  +  \  p"1/p  h.  sup  |b(x)|, 

J  t  xel •  1  J  xel . 

J  J 


1 


h. 


J 


r  j  (x)dx , 
J 


t 


we  now  introduce  the  error  corrector 


Here  C1 ,  are  constants  dependent  on  p  and  the  bounds  (2.2)  which  can  be 
estimated  from  above.  Under  our  smoothness  assumptions  (2.21)  the 
quantities  5.  are  of  order  h(A)  and  hence  the  first  term  of  (2.23)  is  of 
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order  c(a)  h(A).  With  the  aid  of  the  basic  a-priori  relation  between 


Hp(D 


j e 1 1  n  and  because  u'(a)  is  constant  on  I.,  it  can 


Hp<» 


be  shown  that  also  the  second  term  of  (2.23)  is  of  the  order  e( a)  h(A). 

This  leads  to  the  following  result: 

Theorem  2.3:  If  a’,b,f  e  Lp(I)  then  <J>(A)  is  an  error  corrector.  Moreover, 
if  the  conditions  (2.21)  hold,  then 


(2.24) 


!e(A)i 


0  ,  as  h(A)  -*-0,  A  e  P(k)  , 


These  are  only  some  examples  of  possible  error-indicators,  estimators, 
and  correctors  which  have  the  desired  properties.  There  is  no  unique  form 
for  these  quantities  and  various  other  definitions  are  possible. 


2.4  A  Numerical  Example 


We  consider  the  problem 


(2.25) 


-((x+t)1/10  u1)'  +  u  =  f ( x ) 


u(0)  =  u(l )  =  0,  t 


where  f  is  chosen  such  that  the  exact  solution  has  the  form 
(2.26)  uQ(x)  =  (x+t)1/2  -  [(l-x)t1/2  +  x(l+t)1/2],  x  e  I. 


This  allows  us  to  compute  here  the  efficiency  index  (2.14). 


Tables  2.1  and  2.2  below  show  the  errors  in  the  L9-  and  LQ-norm, 

c.  o 

respectively,  for  uniform  partitions  on  I  with  h=l/20,  h=l/40,  and 
h  =  l/80,  as  well  as  for  the  non-uniform  partitions  20A  with  m=20  given 
in  Table  2.3.  The  error  indicators  (2.17)  and  estimators  (2.19)  were 
used.  The  partitions  20A  are  optimal  in  the  sense  defined  in  Section  2.5 
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j 

- ?t - 

j 

j 

_ 

X 

A 

p=2 

P“8 

p*2 

p=8 

0 

o 

0 

10 

.2397 

.1705 

1 

. 1 326 ( - 1  ) 

.91 55 ( -2 ) 

11 

.2828 

.2042 

2 

. 2822 ( - 1 ) 

. 1 942 ( - 1  ) 

12 

.3312 

.2435 

3 

. 4508( - 1 ) 

. 3096 ( - 1 ) 

13 

.3856 

.2895 

4 

. 6408 ( - 1 ) 

. 4398 ( - 1 ) 

14 

.4465 

.3438 

5 

. 8546 ( - 1 ) 

. 587  3 ( - 1  ) 

15 

.5148 

.4083 

6 

.1095 

.  7547 ( - 1  ) 

16 

.5913 

.4854 

7 

.1336 

.9456(-l ) 

17 

.6769 

.5783 

8 

.1670 

.1164 

18 

.7228 

.6910 

9 

.2012 

.1415 

19 

.8800 

.8291 

_ 

20 

1 

1 

Table  2.3:  The  non-uniform  meshes  20A. 

These  results  suggest  the  following  observations: 

(a)  For  uniform  partitions  the  asymptotic  decrease  of  the  error  is  linear 
in  h  under  both  norms. 

(b)  For  the  1,,-norm  the  efficiency  index  o  converges  to  1  with  |je||^ 

c  r  ^ 

t 

(c)  The  efficiency  index  improves  with  decreasing  relative  error  and,  in 
particular,  it  is  better  for  the  non-uniform  optimal  partitions  than  for 
the  corresponding  uniform  ones. 

(d)  The  error  estimators  are  more  effective  under  the  L.,-norm  than  under 


the  Lg-norm. 
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Observations  such  as  these  are  typical  for  much  more  general  problems 
than  the  present  one.  In  particular,  as  we  shall  see  in  the  next  chapter, 
they  also  apply  to  the  two-dimensional  case. 


2.5  Optimal  Meshes 


We  restrict  our  consideration  now  to  a  special  class  of  partitions. 
For  this  let  U  c'(l)  be  any  strictly  increasing  function  with  £(0)  =  0 
and  £(1)  =  1.  Then  a  partition  (2.6)  is  called  a  ( c,ni) -part i tion  and 
denoted  by  AU,m)  if 


(2.27) 


C(xj  )  *  *■  .  j  =0 , 1 . n  =  in (  a  ) 


The  set  of  all  U,in)-partitions  A(£,m)  associated  with  a  fixed  £  and  varying 
m  is  called  the  family  of  graded  partitions  induced  by  the  grading  function 

r 

For  such  families  of  graded  partitions  it  is  possible  to  show  --  under 
suitable  conditions  on  uQ  and  £  --  that 


(2.23) 


U)(p,  Uq,  £)(1  +  0  (  1  )  ) 


This  will  allow  us  to  associate  with  every  member  A(c,m)  of  the  family  the 
asymptotic  error  <o(p,  u^,  £)/m.  Moreover,  it  turns  out  that  there  exists 
a  grading  function  such  that  for  all  other  grading  functions  t  we  have 

(2.1.9)  w(P,Uq,  !-q)2Luj(p*  0  • 

This  means  that  £g  induces  an  asymptotically  optimal  family  of  graded  partitions. 


Our  discussion  here  will  be  restricted  to  the  case  p  =  2.  For  all 
details  we  refer  to  [6  1. 

The  basic  result  may  be  phrased  as  follows: 

Theorem  2.4:  Suppose  that  the  conditions  (2.21)  hold.  Then  there  exists  a 

5 

grading  function  Cg  such  that  m)  e  P(y)  for  all  m  ^  1  and  that  for  any 
other  grading  function  s  f  for  which  A(i.,m)  £  P(<)  for  some  1  <_  «  <  2 
we  have 

(2.30)  |  |e(  a(  Cg  m) )  |  lE>2  <  I  le(A(  f-n’))  I  lE>2  .  vm>m0U). 

Moreover,  if  £  *  +  \4>  is  a  grading  function  for  all  sufficiently  small 

X  then 

(2.31 )  |  |eUU0  +  U  ,  m) )  |  |£  2  <  |  |e(d(  ,  m  ) )  1 1 E  2  +  0(\2) ,  as  X  ■*  ® 

for  all  sufficiently  large  m. 

In  most  practical  problems  the  optimal  partition  is  rather  sensitive 
to  small  changes  of  the  data.  But  (2.31)  shows  that  the  optimal  error  changes 
only  with  the  square  of  any  change  of  the  optimal  partition  and  hence  is 
considerably  more  stable  than  the  latter.  In  other  words,  there  is  no  reason 
why  we  should  attempt  to  compute  the  optimal  partition  at  high  cost  when 
near-optimal  partitions  produce  almost  the  same  error  as  the  optimal  one. 

The  optimal  grading  function  Cg  depends  on  Ug  and  p  and  is  not  difficult 
to  obtain  when  Ug  is  known.  In  fact,  the  partitions  of  Table  2.3  are 
Uq,  20)-partit1ons  for  (2.25)  and  were  computed  from  the  exact  solution 
(2.26).  The  table  clearly  shows  the  dependence  of  ^  on  the  value  of  p. 
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In  general  ,  when  Uq  is  unknown,  it  turns  out  that  we  may  approximate  ^  by 
means  of  the  following  equilibration  principal  for  the  error  indicators: 

Theorem  2.5:  Let  the  conditions  (2.21)  hold  and  suppose  that  for  some  family 
T  of  partitions  (2.6)  we  have 

(2.32)  Hj  ( A)  =  u(l  +  Cl(h(A)p),j  (a),  as  F(a)  •>  0,  A  c  T  , 

with  p  =  1/12.  Then  T  c  P(|)  and 

(2-33)  I (e(A) ( (E  2  <  ((e(A(c0,m))((E>2  (1  +  0(h(A)p),as  F(a)  -  0,  a  e  T. 

From  a  practical  viewpoint  theorems  2.4  and  2.5  suggest  the  use  of 
adaptive  techniques  for  the  simultaneous  construction  of  nearly  optimal 
partitions  and  their  error  estimators.  This  will  be  discussed  in  Chapter  4 
below. 


3. 


Two-Dimensional  Linear  Problems 


3 • 1  The  Model  Problems 

As  before,  for  simplicity  of  discussion  we  present  the  results  for 
certain  model  problems.  More  specifically,  we  consider  the  following  two 
elliptic  boundary  value  problems  on  the  unit  square 


(3.1) 


ft  ■  {x. 


e  R 


0  <  Xj ,  x,,  <  1 ) 


with  Dirichlet  boundary  conditions: 


(a)  Problem  L: 
(3.2) 


.2  , 

AU  +  p  -----  *  0, 

y  dx:)y 


on  ft,  0  p  <  2 


(3.3)  u  *  g  ,  on  3ft 

Here  the  function  g  is  chosen  such  that  the  exact  solution  of  (3.2),  (3.3) 
is  given  by 


1/2 

Uq  =  r  cos  0/2 


(3.4)  r  =  t  2+p  (X1  "  1  +  x2^  +  2-p  ^X1  '  1  “  x2^ 


I  .  /T+n  I  X1  ~  t  “  Xp  A 
0  =  arctan  I  r - »  al~"~  1  *  t  >  1 . 


p  X,  -  t  +  X, 

r  1  c 


Hence  for  p  =  0,  r  and  0  are  polar  coordinates  centered  at  (t,0). 
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(b)  Problem_E: 

2  1 

(3.5)  V  u  +  v  •  vu  =  0  on  ft  ,  u  =  (u^  ,  u9) , 

(3.6)  u»£,on3H,0<io<  ^  ,  £  =  (g1 ,  g2) . 

Here  ^  is  chosen  such  that  the  exact  solution  of  (3.5),  (3.6)  is  as  follows: 
u1  *  j  r^2CU-l )  cos  ^  -  (k-1  )  sin  ^  -  j  cos  |g  -  ^  sin  |-Q  ] 

(3.7) 

U£  =  j  r1/2[(K+l)  cos  £  +  U+l  )  sin  j  +  ^  cos  |e-^sin  ^  0  ] , 

where  <  =  3-4o  and  r,G  are  the  polar  coordinates  centered  at  (t.O),  t  >  1. 

For  t  *■  1  the  solution  has  a  singularity  with  a  character  typical  for  a  corner. 

Both  problems  have  two  parameters,  namely,  t,p  in  the  case  of  problem  L 
and  t,j  in  that  of  problem  E.  Here  t  controls  the  regularity  of  the  solution 
while  p  and  0  characterize  the  degeneracy  of  the  problems. 

We  use  admissible  partitions  A  of  the  unit  square  (3.1)  which  are 
defined  recursively  by  the  following  two  rules: 

(i)  The  undivided  square  is  an  admissible  partition  of  itself. 

(ii)  If  A  is  an  admissible  partition  of  (3.1)  and  s  represents  any 
square  of  A  ,  then  a  new  admissible  partition  A’  is  obtained  if 

(3.8)  s  is  subdivided  into  four  congruent  squares  of  half  the  side- 
length  of  s. 

On  these  admissible  partitions  we  use  the  finite  element  method  with 
test  function  which  are  continuous  on  ft  and  bilinear  on  each  square  of 
A.  Besides  the  usual  nodal  points,  here  called  regular  points,  these 
partitions  contain  nodal  points  inside  a  which  are  not  in  the  corner  of 


all  the  squares  incident  with  them.  These  points  are  called  irregular 
points  and  are  marked  by  an  x  in  Figure  3.1.  Since  the  test  functions  are 
continuous  on  ,  their  values  at  these  irregular  points  are  not  free 
but  are  specified  by  continuity.  For  instance,  the  value  at  point  B  is 
the  average  of  the  values  at  C  and  D. 

As  usual,  with  any  regular  nodal  point,  say  "i",  we  associate  the 
standard  basis  function  and  denote  its  support  by  i :. .  For  instance,  the 
support  associated  with  the  nodal  point  A  of  Figure  3.1  is  the  indicated 
shaded  set. 

For  any  admissible  partition  A  let  the  irregularity  index  x(a)  be 
defined  as  the  maximal  number  of  irregular  points  between  any  two  adjacent 
regular  points.  For  instance,  for  the  partition  of  Figure  3.1  we  have 
X(A)  =  4.  We  shall  always  restrict  our  consideration  to  a  family  P  of 
admissible  partitions  with  \{a)  <_  \  <  *>  for  all  A  e  P. 

3.2  Error  Indicators  and  Estimators 

For  two-dimensional  problems  the  theory  of  a-posteriori  estimates  is 
not  yet  as  fully  developed  as  for  the  one-dimensional  case.  We  shall  not 
go  into  the  details  of  the  theory  here  but  indicate  only  some  of  the 
typical  theorems  and  present  certain  representati ve  experimental  results. 
As  in  Chapter  2  we  shall  use  the  L^-stress-energy  norms. 

The  basic  result  may  be  broadly  phrased  as  follows  (see  eg  [4 1 ) : 

.  / 

Theorem  3.1:  Let  A  e  P  be  an  admissible  partition  on  Q  and  u ( A )  the 
corresponding  finite  element  solution.  For  each  regular  nodal  point  i  of 
A  define  e  as  the  exact  solution  of  the  given  problem  on  the 


2  0 


associated  support  set  0..  which  satisfies  the  boundary  conditions 
z ^  =  u(a)  on  an.  and  =  g  on  an^  n  an. 

Moreover,  set 


(3-9)  ni  Mzi  I  I E.2,n1  ’ 

Then  there  exist  constants  >_  K1  >0  such  that 


(3.10) 


where  e  =  uQ  -  u(a)  is  the  actual  error.  The  constants  ,K^  depend  on  the 

problem  and  P,  but  not  on  A  and  u(a). 

This  shows  that  once  again  an  upper-  and  lower-estimator  can  be  obtained 
on  the  basis  of  strictly  local  information.  Of  course,  in  practice  the  n. 

are  not  computable.  However,  by  using  appropriate  bounds  it  is  possible  to 
construct  computable  error  indicators  n  •  for  the  elements  of  A  which 

J 

provide  estimators  that  are  asymptotically  equivalent  with  the  estimate  (3.10). 


We  shall  not  enter  into  the  details  of  the  construction  of  these  error 
indicators.  Briefly,  if  s.  is  one  of  the  elements  of  the  given  partition 
(see  Figure  3.2),  then  the  error  indicator  associated  with  s.  is  composed 
of  h-  and  of  the  values  of  u(a)  and  the  jumps  of  the  normal  derivatives  of 
u(a)  at  the  corners  of  s^.  The  particular  form  of  n.  depends  only  on  the 
differential  equation  and  the  specified  norm. 


2 


Figure  3.2 


The  basic  asymptotic  property  of  these  indicators  can  be  seen  from  the 
following  special  result: 

2 

Theorem  3.2:  Let  the  domain  be  all  of  R  .  Then,  for  a  family  of  uniform 
partitions  with  side-length  h(A)  and  for  sufficiently  smooth  exact  solutions 
Uq,  the  estimator 


(3.11) 


satisfies 


:(A)  =  [  V  n?  ]1/2 


(3.12)  I  lu0  -  u( a)  1 1^2  =  e(A)2  +  0(h(A)4),  as  h(a)  >  0. 

In  all  reasonable  cases  we  expect  that  llell^  2  -  ^(a)  »  whence 


|U0  -  u(A)ME>2  =  c( A)(l  ♦  0(h(A)^)), 


=  e(A)  (1  +  0(  e  ( A)4" ) ) ,  as  h(A)  -  0. 


(3.13) 


Hence  the  efficiency  index  (2.14)  will  converge  to  one  with  the  square 
of  the  error. 

Nearby  optimal  meshes  may  now  be  constructed  adaptively  as  in  the 
one-dimensional  case  by  means  of  the  equilibration-principal  for  the  error- 
indicators.  Because  a  number  of  general  mathematical  questions  are  still 
open,  we  will  not  go  here  into  further  theoretical  details.  Instead  we 
present  various  experimental  results  which  support  a  number  of  conjectures 
for  which  there  is  as  yet  no  sufficiently  general  proof. 

Conjecture  1 :  For  sufficiently  smooth  solutions  on  a  bounded  domain  and  for 
uniform  partitions,  the  efficiency  index  with  respect  to  the  L0-stress 
energy  norm  converges  to  one  with  the  square  of  the  error. 

This  is  clearly  seen  in  Figure  3.3  for  the  solution  of  problem  E  in 
the  two  cases 

case  E, :  t  =  2 . 0  o  =  0.0,  uniform  partition, 

(3.14)  1 

case  E^:  t*2.0,o  =  0.45,  uniform  partition. 

The  figure  shows  in  log/log-scale  the  square  of  the  error  (under  the 
U -stress  energy  norm)  as  well  as  the  deviation  from  one  of  the  efficiency 
index  0  as  a  function  of  the  number  N°  of  degrees  of  freedom  of  the  finite 
element  solution. 

Conjecture  2:  For  uniform  partitions  on  bounded  domains  the  rate  of 
convergence  of  the  efficiency  index  to  one  is  more  strongly  influenced  by 
the  smoothness  properties  of  the  exact  solution  uQ  than  by  the  convergence 


of  the  finite  element  solutions  to  u^. 
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This  is  supported  by  Figure  3.4.  It  has  the  same  form  as  Figure  3.3 
and  shows  results  for  problem  E,  case  E^  and 

(3.15)  case  E^:  t  =  1.05,  o  =  0.0,  uniform  partition, 

as  well  as  for  problem  L  in  the  two  cases 

case  L-j :  t  =  1.05,  p  =  0.0,  uniform  partition, 

(3.16) 

case  L9:  t  =  1.001,  p  =  0.0,  uniform  partition. 

Evidently,  the  singularity  becomes  stronger  as  t  tends  to  one. 

The  next,  very  important  conjecture  concerns  the  adaptively  constructed 
partitions.  Linder  suitable  conditions  it  is  known  that  for  non-smooth 
solutions  there  exist  sequences  of  properly  refined  meshes  for  which  the 
rate  of  convergence  of  the  error  is  comparable  to  that  for  smooth  solutions, 
provided,  of  course,  the  same  type  of  elements  are  used.  In  this  light  the 
conjecture  is  not  surprising: 

Conjecture  3:  For  the  adaptfvely  constructed  sequences  of  meshes  the 
rate  of  convergence  of  the  error  (under  the  L9-stress  energy  norm)  and  of 
the  efficiency  index  compares  with  that  for  smooth  solutions  and  uniform 
partitions.  In  particular,  the  efficiency  index  converges  essentially 
with  the  square  of  the  error. 

Figure  3.5  gives  results  for  problem  L  in  the  two  cases 

case  L9:  t  =  1.05,  p  =  0.0,  adaptively  refined  partition, 

(3.17) 

case  Lj:  t  =  1.05,  p  =  0.0,  adaptively  refined  partition. 


MP* 
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A  comparison  of  Figure  3.4  and  3.5  shows  the  contrast  between  the  solution 
of  the  same  problem  using  uniform-and  adaptively-refined  partitions.  The 
behavior  of  the  efficiency  index  close  to  one  is  influenced  by  various 
factors  and  is  not  fully  understood.  Some  aspect  of  this  is  the  topic  of 
the  next  conjecture. 
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Conjecture  4:  For  completely  arbitrary  partitions  of  fixed  irregularity 
index  the  efficiency  index  does  not  converge  to  one  but  has  reasonable 
bounds >  dependent  on  the  regularity  index. 

The  partitions  of  Figure  3.6  were  randomly  generated.  The  conjecture 
is  certainly  supported  by  the  results  of  Table  3.1  for  problem  L  in  the 
case 

(3.18)  case  L^:  t  =  2.0,  p  =  0.0,  random  partitions. 


Partition 

No  of 
nodes 

No  of 
elements 

0 

1 

9 

16 

.976 

2 

21 

40 

.868 

3 

34 

64 

.871 

4 

49 

88 

.868 

5 

72 

115 

.926 

6 

103 

166 

.901 

7 

121 

190 

.894 

8 

132 

214 

.877 

Table  3.1  Efficiency  indices  for  random  partitions 
As  before  the  efficiency  index  0  was  computed  under  the  1,,-stress-energy  norm. 
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Conjecture  5:  The  values  of  the  parameters  p  and  a  controlling  the  degeneracy 
of  the  solutions  of  (3.2),  (3.3)  and  (3.5),  (3.6),  respectively,  do  not 
significantly  influence  the  efficiency  index. 

Conjecture  6:  The  adaptive  construction  of  the  partitions  is  effective  both 
for  the  1*2  and  the  L^-energy  norms.  The  efficiency  index  for  the  Lg-norm  is 
better  than  that  for  the  L  -norm,  but  in  both  cases  it  is  acceptable. 

oo 

For  problem  L  and  the  case  of  (3.17),  Figures  3.7  and  3.8  show 
the  adaptively  constructed  meshes  for  the  l^-energy  norm  and  the  L^-energy  norm, 
respectively.  Tables  3.2  and  3.3  summarize  the  errors  and  efficiency  indices 
for  these  cases.  Here  the  L^-energy  norm  was  computed  as  the  maximum  value 
over  the  set  of  the  four  Gaussian  points  in  each  element.  The  table  headings 
are  as  follows: 

(1)  Partition  identification  (5)  Square  of  the  exact  error 

(2)  Number  of  degrees  of  freedom  (6)  Relative  error  in  percent 

(3)  Number  of  elements  (7)  Efficiency  index 

(4)  Square  of  the  error  estimator 
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4.  Adaptive  Solution  of  Linear  Problems 


4.1  The  Basic  Design  of  FEARS. 

In  recent  years,  many,  often  large,  general-  and  special  purpose 
software  systems  for  practical  finite  element  computations  have  been 
developed  and  are  used  extensively.  The  same  time  has  brought  considerable 
advances  in  the  theoretical  analysis  of  the  method,  the  development  of 
efficient  solution  procedures,  and  the  understanding  of  modern  software 
systems.  There  appears  to  be  a  need  for  utilizing  these  and  related 
advances  in  the  design  of  the  next  generation  of  finite  element  software. 

Currently  an  experimental  software  system  is  under  development  at  the 
University  of  Maryland  at  College  Park,  MD.  It  has  been  named  FEARS 
(F_inite  Element  Adaptive  Research  Solver),  and  has  the  following  principal 
design  properties. 

(i)  The  system  constitutes  an  appl ications-independent  finite 
element  solver  for  a  particular  class  of  two-dimensional, 
linear  elliptic  boundary  value  problem. 

(ii)  The  problems  are  defined  by  a  weak  mathematical  formulation  in 
terms  of  a  bilinear  form. 

(iii)  The  a-posteriori  error  estimates  discussed  in  Chapter  2  above 
are  used  to  control  the  reliability  of  the  solution.  Different 
norms  may  be  chosen  by  the  user. 

(iv)  Adaptive  mesh-refinement,  based  on  the  error  estimates,  are  to 
provide  for  near  optimal  meshes  within  a  prescribed  cost  range. 

(v)  The  systems-design  takes  advantage  of  some  of  the  natural 
parallelism  and  modularity  of  the  finite  element  method. 


The  system  represents  an  experimental  prototype  rather  than  a  production 
system.  Accordingly,  extensive  provisions  for  performance  evaluation  are 
incorporated  in  it.  A  general  description  of  the  overall  design  is  given  in 
[7].  A  first,  preliminary,  and,  as  yet,  limited  version  of  the  system  has 
been  operational  for  about  a  year. 

As  provided  by  the  second  design  property,  the  system  is  based  on  a  weak 
formulation  of  the  given  problem  in  terms  of  an  appropriate  bilinear  form  B 
on  certain  Hilbert  spaces  ,  H?  and  a  corresponding  (load)  functional  f  on  H9 
In  other  words,  the  underlying  mathematical  problem  is  to  obtain  the  (unique) 
Uq  e  H.j  such  that 

(4.1)  B(u0,v)  =  f(v),  \/v  e  H?. 

The  finite  element  method  then  proceeds  to  the  construction  of  suitable  finite 
dimensional  subspaces  M-|  c  ,  M0  c  and  the  numerical  computation  of  the 
solution  u(a)  c  M.j  of 

(4.2)  B(u(A),  v)  =  f  ( v ) ,  vv  e 

The  choice  of  B  and  f  determines  the  problem  and  the  structure  of  the  finite 
element  method.  The  selection  of  the  spaces  ,  H?  affects  the  norms  used  in 
the  error  analysis  and  hence  the  adaptive  control  of  the  overall  solution 
process.  Finally,  the  spaces  derive  from  the  construction  of  the 

particular  meshes  and  the  selection  of  the  specific  elements.  We  refer  to 
[  8]  for  the  theoretical  background  of  this  formulation  and  to  [  2 1  for 
some  of  the  reasoning  behind  its  use  as  a  basis  for  the  systems-design. 

The  specific  form  of  B  and  f  used  in  FEARS  may  be  found  in  [7].  We 
note  here  only  that  the  selected  formulation  is  fairly  general.  It  includes 
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not  only  most  classical  problems  and  finite  element  methods  but  also  the 
elliptic  problems  of  Dougl i s-Ni renberg  type,  the  forms  arising  in  the  use 
of  the  least-squares  method  for  first-order  systems,  various  versions  of  the 
Lagrange  multiplier  methods,  and  many  others.  Suitable  pre-processing 
facilities  can  be  introduced  to  simplify  the  problem  definition. 

A  natural  parallel  process  structure  for  the  system  derives  from  the 

familiar  substructure  analysis  in  engineering  design.  For  this  purpose, 

_  ? 

the  domain  a  c  R  is  defined  as  the  union 

(4.3)  Q  *  Oj  U  n2  U...U  ^ 

_  2 

of  finitely  many  compact  subsets  P..j  c  R  which  have  non-empty  interiors 

il.  that  are  disjoint,  il.  H  n.  =  0,  i  /  j .  Each  subdomain  sT,.  is  assumed  to 

2 

be  the  diffeomorphic  image  of  some  fundamental  figure  F  in  R  ,  (see  Figure  4.1). 
On  the  intersections  of  the  subdomains  these  di ffeomorphims  *.:F  -*■  a. 

have  to  satisfy  appropriate  compatibility  conditions;  for  details  we  refer  to 
[  7  ]• 

The  finite  element  meshes  on  each  n.  consist  of  curvilinear  elements 
which  are  first  defined  on  F  and  then  mapped  by  <!>..  into  fi. .  Thus  the 
mesh-construction  takes  place  in  F.  In  order  to  support  the  adaptive  mesh- 
refinement  process,  the  fundamental  figure  should  admit  the  definition  of 
a  simple  hierarchy  of  partitions  consisting  of  geometrical  figures  of  the 
same  shape,  and  there  should  be  a  simple  subdivision  algorithm  which  allows 
for  the  construction  of  the  partitions  of  the  hierarchy. 

For  the  design  of  FEARS  the  fundamental  figure  F  was  chosen  to  be  the 
closure  Qq  of  the  open  unit  square 

Q0  =  {x  e  R2  0  <  x.  <  1 ,  1-1,2} 


(4.4) 
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As  in  the  example  of  Section  3.1,  the  admissible  partitions  A  on  Qq  are  then 
defined  recursively  by  the  two  rules  (3.8).  The  current  version  of  FEARS 
incorporates  only  bilinear  conforming  elements  on  the  squares  of  these 
partitions  of  F  which  are  then  mapped  into  the  subdomains  ft.. .  This  restric¬ 
tion  of  the  choice  of  elements  was  dictated  only  by  limitations  of  available 
programming  support. 

4.2  Computational  Aspects 

The  introduction  of  parallelism  indicated  in  the  design  property  (v) 
appears  to  be  a  particularly  novel  aspect  of  the  design  of  the  FEARS  system. 
This  parallelism  is  on  the  procedural  level  rather  than  the  instruction- 
level;  that  is,  it  is  defined  in  terms  of  processes  which  are  autonomous 
units  with  their  own  programs  and  data.  These  processes  run  in  parallel 
and  communicate  asynchronously  in  a  limited  and  highly  structured  manner. 

A  basic  design  decision  was  not  to  store  data  redundantly  whenever  possible. 
Thus  each  process  contains  almost  all  the  information  needed  for  it.  The 
use  of  parallel  processes  in  the  design  makes  it  possible  to  apply  multiple 
processors  effectively.  But  this  does  not  limit  the  system's  use  to  such 
specialized  architectures;  the  operating  system  of  a  typical  large  computer 
simulates  parallel  processes,  and  the  data  segmentation  should  provide 
significant  advantages. 

There  are  four  classes  of  processes  in  FEARS,  namely,  the 

(i)  geometrical  processes 

(ii)  user  processes, 

(4.5) 

(iii)  linear  systems  solvers,  and  the 

(iv)  control  process. 
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Since  the  closed  subdomains  ft.,  of  (4.3)  intersect  and  hence  introduce 
redundant  data,  we  are  led  to  work  exclusively  with  open  subdomains. 
Accordingly,  a  separate  process  is  associated  with  each  one  of  the  2-subdomains 
ft-|,...,  ft^,  with  each  one  of  the  relatively  open,  one-dimensional  intersections 
(called  1 -subdomains)  of  two  of  the  ft.,  and  with  each  O-dimensional  inter¬ 
sections  (called  O-subdomains)  of  several  of  them.  The  geometric  processes 
contain  all  the  information  relevant  to  that  particular  subdomain.  This 
includes  the  bilinear  form,  load  functional,  mapping  between  Qq  and  the 
subdomain,  boundary  conditions  where  they  apply,  as  well  as  the  current  mesh 
data  and  associated  solution.  The  only  data  set  stored  more  than  once  is  the 
adjacency  relation  between  the  subdomains. 

A  problem  is  defined  by  a  given  domain  (4.3),  bilinear  form  B,  and 
associated  error  norms.  Each  instance  of  the  system  is  initialized  with 
a  particular  problem.  It  can  accommodate  several  users  working  simultaneously 
on  this  problem,  each  one  with  his  own  boundary  conditions,  load  functional, 
sequence  of  meshes,  and  solutions.  The  user  processes  contain  the  data 
particular  to  that  user  and  any  information  needed  for  post-processing, 
graphics,  etc. 

The  third  class  of  processes  represent  the  linear  systems  solvers  which 
have  no  permanent  data.  Finally,  the  efforts  of  the  various  processes  are 
orchestrated  by  a  central  control  process.  It  receives  commands  from  the 
current  user,  gives  commands  to  the  geometric  processes  and  the  system 
solvers,  receives  reports  back  from  them,  and  sends  reports  back  to  the  user. 

Basically,  after  some  initialization  phase,  a  user  gives  a  command  which 
tells  the  system  to  refine  the  current  mesh,  and  to  compute  the  corresponding 
finite  element  solution  so  as  to  achieve  a  prescribed  error  tolerance  without 
exceeding  a  given  computational  cost. 


Hence,  if  we  disregard  parallel ity  for  the  moment,  a  user  requests  the 
system  to  execute  the  following  algorithm: 

1.  [Input]  'initial  mesh'  a;  'tolerances'  <5  kS,  'Vel’ 

'maximal  cost'  y  ; 

max 

2.  [Initialization]  'cost'  y :*0;  'cost  increment'  Ay:=0; 

'error  estimator'  >.  'norm'  v(a):s>»; 

'cut-off  parameter'  n:"0; 

3.  [Main  Loopl 

While  (y  +  Ay  •-  y  )  and  (.  (a)  >  A  .  +  <$  ,  v(a)  )  do 

'  '  'max  '  abs  rel  — 

3.1  assemble  stiffness  equations;  compute  solution  u(a); 
compute  norm  v(a)  *  ||u(a)||;  update  cost  y; 

3.2  for  each  element  t  .  A  compute  error  indicator  n(r); 
accumulate  estimator  .(a);  if  A  is  not  the  initial  mesh 
determine  the  cut-off  parameter  p ; 

3.3  subdivide  each  element  r  .  a  with  n(t)  >  n;  for  the 
new  elements  r'  obtained  from  v  set  r^(r')  3  n(r); 

3.4  predict  cost  Increment  Ay  of  next  solution  phase. 

briefly,  the  algorithm  continues  until  either  the  maximum  allotted  cost 
will  be  exceeded  or  the  error  estimator  is  sufficiently  small.  Except,  for  the 
initial  in^sh  the  cut-off  parameter  n  is  the  maximum  value  of  the  quotient 
n(  r )  /rig( T )  where  hq(i  )  is  the  error  indicator  of  the  element  from  which  i 
was  generated  by  subdivision.  This  represents  a  simple  procedure  to  force  the 
error  indicators  to  be  closer  together. 
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In  the  procedure  we  should  distinguish  between  the  mesh  design  phase 
and  the  computation  of  the  final  solution.  During  the  construction  of  the 
various  intermediate  partitions  A  there  is  no  need  to  spend  too  much  effort 
on  the  solutions.  It  turns  out  that  the  sequence  of  refined  meshes  is  not 
very  much  affected  if  the  new  solution  values  are  only  computed  on  the  newly 
subdivided  portions  of  the  mesh;  that  is,  if  on  the  other  parts  of  the  mesh 
the  old  solution  is  retained  and,  where  needed,  used  as  boundary  condition. 

We  even  found  that  for  two  to  three  passes  it  is  possible  to  compute  the  effect 
of  the  refinement  on  each  element  separately.  Thus  in  FEARS  a  combination  of 
these  types  of  "short"  solution  passes  is  used  consistently  except  on  the  first 
two  meshes,  where  complete  solutions  are  needed  to  start  the  procedure,  and 
on  the  final  mesh. 

So  far  we  disregard  the  effects  of  parallelism.  Generally,  the  following 
parallel  activities  take  place: 

(a)  Several  users  perform  pre-  and  post-processing  while  one 
user  is  in  the  solution  phase. 

(b)  The  2-subdomain  processes  compute  their  error  indicators  in 
parallel . 

(c)  The  2-subdomain  processes  refine  their  meshes  in  parallel. 

(d)  The  2-subdoniain  processes  set-up  their  blocks  of  the  stiffness 
matrix  in  parallel. 

(e)  The  2-subdomain  processes  decompose  their  blocks  of  the 
stiffness  matrix  in  parallel  and  send  updates  to  a  system 
solver  process. 

(f)  During  short  solution  passes,  "disjoint"  solution  sets  are 
computed  in  parallel. 

In  connection  with  the  activities  (d)  and  (e)  above  it  should  be  noted 
that  the  macro-stiffness  matrix  has  the  form 


■ 
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B 
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where  each  diagonal  block  corresponds  to  a  2-subdomain  and  the  blocks 

C^,  B  of  the  border  contain  the  data  for  the  0-  and  1 -subdomains . 

Clearly,  each  2-subdomain  process  may  perform  the  following  tasks  in 
parallel : 

(i)  Generate  the  blocks  Ai ,  Ci  and  send  the  corresponding 

contribution  for  the  block  B  to  a  separate  solver  process. 

(ii)  Decompose  A.,  modify  accordingly,  and  send  the  resulting 
modification  of  B  to  the  solver. 

After  all  2-subdoamin  processes  finish  these  steps,  the  solver  process  completes 
the  decomposition  of  B  and  now  the  new  solution  is  obtained  in  an  obvious  way. 
Some  modi fications ,  of  course,  are  needed  in  the  case  of  short  solution  passes. 

The  efficiency  of  the  refinement  strategy  depends  critically  on  the  design 
of  the  data  structure  for  the  meshes.  A  widely  used  data  structure  for  finite 
element  computations  is  based  on  a  list  of  the  nodes  each  pointing  to  the 
elements  to  which  it  belongs,  and  a  list  of  the  elements  which  in  turn  point 
to  the  nodes  incident  with  them.  This  essentially  static  structure  is  not  very 
efficient  when  mesh  refinements  are  introduced.  Instead  the  recursive  defini¬ 
tion  (3.8)  suggests  the  use  of  a  tree  structure  that  corresponds  to  the 
refinement  process.  The  tree  structure  developed  for  and  used  in  FEARS  is 
described  in  [  9]  and  we  sketch  here  only  the  basic  idea. 

As  discussed  before,  the  mesh  construction  takes  place  on  the  closure 
of  the  square  (4.4).  Hence  the  tree  structure  represents  the  refinement 


JL. 
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process  on  Qq  as  defined  by  (3.8).  As  mentioned  in  Section  3.1,  in  the 
resulting  partitions  we  need  to  distinguish  between  regular  and  irregular 
nodal  points  (see  Figure  3.1).  The  nodes  of  the  tree  correspond  to  all 
regular  nodal  points  and  to  the  center  of  all  undivided  squares  of  the 
partition.  The  successor  relation  on  the  tree  represents  the  subdivision 
process.  For  details  we  refer  to  [ 9  1 .  An  example  of  the  resulting  tree 
structure  is  given  in  Figure  4.2. 

The  nodes  of  the  tree  carry  labels  of  the  form  (X-j,  \0), 

X-| ,  \ 2  =  +1  -  which  permit  the  retrieval  of  all  needed  geometrical 

information.  In  fact,  as  shown  in  [9],  simple  algorithms  are  available 
to  perform  the  following  tasks: 

(a)  For  a  given  tree-node  representing  a  divided  or  undivided  square 
s  determine  --  as  far  as  they  exists  -- 

(al )  the  corners  of  s 

(a2)  the  midpoints  of  the  sides  of  s 

(a3)  the  neighboring  squares 

(b)  For  a  terminal  tree  node  representing  an  undivided  squares 

(bl)  subdivide  s  and  extend  the  tree 

(b2)  compute  the  elementary  stiffness  matrix  of  s 

(c)  Decompose  the  block  A.,  of  (4.6)  defined  by  the  tree 

using  a  natural  point  order  induced  by  the  tree  corresponding 
to  nested  dissection.  Thus  there  is  no  need  for  any  renumbering 
algorithms  or  special  strategies. 

In  connection  with  (b2)  we  should  note  that  a  particular  square  s  may 
have  up  to  three  corners  that  are  irregular  nodal  points  of  the  partition. 


2 


6 
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This  has  to  be  taken  into  account  in  the  computation  of  the  corresponding 
elementary  stiffness  matrix.  The  procedure  is  particularly  simple  for  our 
case  of  bilinear  elements.  Let  M  denote  the  standard  elementary  stiffness 
matrix  of  s  when  all  corners  are  regular.  If  some  of  the  corners  are 
irregular,  there  exists  a  --  possibly  rectangular  --  interpolation-matrix 
P  such  that  the  actual  elementary  stiffness  matrix  of  s  has  the  form  P^MP. 
The  matrix  P  depends  only  on  the  geometry  of  the  partition  and  not  on  the 
finite  element  problem  itself.  Thus  P  can  be  computed  directly  from  the 
tree. 


4.3  A  Numerical  Example 

Extensive  numerical  experiments  have  been  conducted  with  the  pilot 
version  of  FEARS.  Here  we  present  only  one  example  involving  the 
solution  of  a  plane,  linear  elasticity  problem  for  an  isotropic, 
homogeneous  material  with  Poisson  ratios  o  =  0  and  a  =  0.45.  The  domain 
is  the  L-shaped  set 


si  =  { x  £  R 


1  2 

0  <  X|  <_  1  ,  0  ^  x^  U  ix  c  R 


°i’i4 


x2  <  11, 


The  displacements  in  the  x^  and  x0  directions  are  denoted  by  u^  and  u^, 
respectively,  and  the  following  boundary  conditions  are  used 

1  3 

u1  =  0,  u^  =  +1  for  4  i  X]  £  4  »  =  0 

u-j  =  0,  u^  =  -1  for  0  <  <  ^,  x0  =  1 . 


The  other  parts  of  the  boundary  ,  are  free  of  stresses.  Thus 
the  problem  may  be  view  as  a  two-dimensional  model  of  an  asymmetric 
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stamp.  We  note  that  the  singularity  of  the  solution  at  (1/4,  0),  (3/4,  0) 
introduces  oscillatory  behavior  in  the  stresses. 

The  refinement  process  begins  with  a  coarse  regular  mesh  with  step-size 
h  =  1/4  (mesh(0))  which  is  then  subdivided  into  a  regular  mesh  of  step  size 
h  =  1/8  (mesh(l)).  For  the  case  o  =  0  the  meshes  (2)  to  (6)  generated  thereafter 
by  the  algorithms  are  shown  in  Figure  4.3  and  4.4.  Moreover,  Figure  4.5  shows 
the  mesh  (7)  for  o  =  0  as  well  as  the  mesh  (5)  for  the  case  a  =  0.45  which 
corresponds  to  it  in  the  size  of  the  error  estimator.  (For  space  reasons  the 
meshes  (2),  (3),  (4)  for  the  case  o  =  0.45  are  not  included).  It  is  interesting 
to  note  the  differences  in  the  effect  of  the  singularities  for  different 
Poisson  ratios  at  these  error  tolerances.  In  all  cases  the  L^-energy  norm  was 
used. 


Tables  4.1  and  4.2  show  some  of  the  results  of  these  computations 


mesh 

ident. 

no  of 
elements 

e(A) 

1  00tt-,4^ — 

1  |U(A) | |E>2 

1 

48 

1 . 073 ( - 1 ) 

10.8':. 

2 

54 

9.890(-2) 

10.1% 

3 

66 

9 . 1 42 ( -2 ) 

9.3% 

4 

78 

8.389(-2) 

8.6% 

5 

84 

7. 977 (-2) 

8.2% 

6 

111 

7 . 233 ( -2 ) 

7.4% 

7 

171 

6.071 (-2) 

6.3% 

Table  4.1 :  "Stamp"  problem,  Poisson  ratio  o  =  0. 


-a 


mesh 

no  of 

eU) 

100riu(i)l iE>2 

ident 

elements 

1 

48 

l.ll(-l) 

25.2% 

2 

84 

9 . 43 ( -2 ) 

21 .8% 

3 

129 

8 . 22 ( -2 

19.3% 

4 

183 

7 . 29 ( -2 ) 

17.2% 

5 

_ 

255 

_ 

6.60(-2) 

15.7% 

Table  4.2  "Stamp"  problem,  Poisson  ratio  o  = 0.45 


Clearly  the  relative  errors  are  considerable  higher  for  o  =  0.45. 
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5.  Non! inear  Problems 

5.1  Error  Indicators  and  Estimators 

The  theory  of  a-postqriori  error  estimates  sketched  in  Chapters  2  and 
3  can  also  be  extended  to  nonlinear  problems.  In  analogy  to  (2.1)  we 
consider  the  elliptic  boundary  value  problem 


-^a(^)  +  b(u)  =  f(x)’  x  e  1  =  C0J1  ’ 

(5.1) 

u(0)  =  u(l)  =  0. 


Clearly  a,b,f  need  to  satisfy  appropriate  conditions,  but  we  shall  not 
enter  into  any  details  here. 

As  before  different  norms  may  be  considered.  For  simplicity  we  shall 
use  the  semi  norm  | u | ^  =  | ] u ' | | L  which  is  here  an  equivalent  norm  on 

o,  ,p  p^ 

HJ(0. 

O 

Let  A  again  denote  the  partitions  (2.6)  of  I  and  S  the  set  of 
continuous  functions  on  I  which  are  linear  on  each  subinterval 
1^  ,  j  =  1 .... ,  m(A),  and  zero  at  the  endpoints  xg  =  0,  =  1.  Then  the 

^  °A 

finite  element  solution  u(a)  e  S  of  (5.1)  is  defined  by 

rl  fl  o. 

(5.2)  j  [a(u')v'  +  b(u)v]dx  =  |  f(x)vdx  ,  V v  e  S'  . 

0  '0 

As  in  the  linear  case,  error  indicators  of  various  form  may  be  defined 

A  .  .A 

on  the  subintervals  I..  We  use  here  the  finite  element  solution  on  I. 

J  J 

when  a  quadratic  term  has  been  added  to  the  computed  solution  u(a).  In 
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other  words,  set  y-  =  u(a)(x.),  j  =  0,...,  m(A)  and 
J  J 


(5.3)  w(x,z)  =  y^_1  +  h  (yj 
Then,  with  the  solution  z  e  r'  of 


-  y.^ )  +  z(x  -  xj_i ) (xj  -  x),  vx 


(5.4) 


Xj 


Xj-1 


a(wx(x,z))(xj._1  +  Xj  -  2x)dx 


fX, 


J  (b(w(x,z))  -  f(x))(x  -  x.  ,)(x- - x)dx=  0, 

J  *  J 


j-1 


we  define  the  error  indicators  by 


(5.5) 


n j  =  1  z j  | »  3  ~  1  >  •  •  •  *  ro( M  » 


and  the  error  estimator  by 


(5.6) 


:(A)  = 


[m{{]  (n$)>] 


1/p 


,  1  <  p  < 


j=l 


max  Hj 


Not  unexpectedly,  (5.4)  represents  a  nonlinear  equation  for  the  error 
indicators.  In  general,  it  is  easily  solvable  by  standard  iterative  processes 
starting  with  the  initial  approximation  z  =  0. 

The  theory  of  these  indicators  and  estimators  is  even  less  developed 
than,  say,  in  the  case  of  the  linear  problems  of  Chapter  3.  But  under 
suitable  smoothness  on  a,b,f  some  of  the  results  of  Chapter  2  can  indeed 
by  carried  over  to  this  nonlinear  case.  This  will  be  discussed  elsewhere. 

. ..  ■-  . . . .  . -  ■  — ■  ■  - .  . . .....  — — ^ 
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5.2  Continuation  Methods 

o 

The  equation  (5.2)  determining  the  finite  element  solution  u(a)  e  S 
of  (5.1)  represents  a  nonlinear  system  of  n  =  m(A)-l  equations  in  the  n 
unknown  y,  =  u(a)(x.),  j  =  l,...,  n.  In  structural  analysis  the  most  common 

J  J 

technique  for  solving  such  systems  is  the  continuation  method,  (see  eg.  [10]). 
For  this  we  assume,  for  example,  that  the  right  side  of  (5.1)  depends  on  a 
parameter  x,  say,  f(x)  =  Xfg(x)  with  a  given  function  fg.  Then  the  resulting 
system  of  equations  has  the  form 

(5.7)  Gy  =  xg,  ye  Rn 

where  G:Rn  Rn  is  given  nonlinear  function  and  g  e  Rn  a  fixed  vector. 

This  system  (5.7)  defines  a  family  of  curves  y  =  y(x)  in  Rn,  and  the  standard 
continuation  approaches  numerically  trace  the  solution  curve  defined  by  the 
initial  condition  y(Xg)  =  y^  where  Gy°  =  xQg.  For  example,  if  a(0)  =  b(0)  =  0 
then  we  may  take  Xg  =  0,  y° =  0. 

It  is  well-known  that  these  methods  encounter  difficulties  when  the 
Jacobain  G'(y)  becomes  singular.  We  have  to  distinguish  here  between  the 
potential  bifurcation  points  where  rank  (G'(y),-q)  <  n  and  the  limit  points 
where  rank  (G'(y),-g)  =  n  >  rank  G'(y).  The  limit  points  do  not  represent 
intrinsic  singularities  of  our  problem.  In  fact,  they  loose  much  of  their 
significance  if  the  solutions  of  (5.7)  are  considered  as  curves  in  Rn+^ . 

For  this  we  add  x  as  (n+l)st  component  to  the  vector  y  and  hence  consider, 
in  general,  an  "underdetermined"  equation 


(5.8) 


Fy  =  c 


involving  a  function  F:Rn*'  »  Rn.  A  regular  C*-solution  of  (5.8)  is  then 
^1  | 

any  curve  In  l<  represented  by  a  C  -parametrizatlon 

(5.9)  y  ■  y(s) ,  scJ  c  R1,  y'  (s)  /  0,  Vs  ,  J 

such  that  F(y(s))  *  c  for  all  s  in  the  open  interval  J.  For  arbitrary  vectors 
c  in  the  range  of  F  the  resulting  collection  of  curves  defines  the  solution 
field  of  F. 

t 

An  analysis  of  such  solution  fields  was  given  in  [11]  and  was  used  there 
to  show  that  a  particular  continuation  method  is  well  defined.  We  shall  not 
enter  into  details  here  but  sketch  only  the  essential  ideas  of  the  approach. 

For  each  vector  y  in  the  regularity  set 

(5.10)  K(F)  -  l.y  .  Rn+1  rank  F'(y)  -  n) 

there  exists  an  unique  vector  z  .  Rnt*  such  that 

(F'(y)\ 

(5.11)  F'(y)z  -  0,  ||z||;,  ■  I.  det  y  T  J  >0 

Hence  the  tangent  mapping 

(5.12)  T:R(F)  *  Rn+1 ,  Tx  -  z 

is  well-defined.  In  [11]  the  following  result  was  proved: 

Theorem  5.1:  Let  F:Rn*^  ►  Rn  be  twice  continuously  differentiable.  Then 
the  tangent  mapping  T  of  (5.12)  is  Lipschitzlan  on  any  compact  subset  of  k’(F), 
and  for  any  y°  ,  0(F)  and  c  *  Fyl1  there  exists  a  unique,  simple,  regular 
(^-solution  in  k’(F)  of  Fy  c  which  has  no  endpoint  in  R(F). 
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These  c' -sol utions  in  R(F)  may  be  parametrized  by  the  solutions  (5.9) 
of  the  autonomous  initial  value  problem 

(5.13)  y*  *  Ty  ,  y(0)  =  y°, 

in  which  case  the  parameter  is  the  arclength. 

The  vector  Ty  represents  the  tangent  of  the  solution  field  at  the  point 
y  e  R(F)  and  y  is  a  limit  point  with  respect  to  the  ith  variable  if  (Ty ) ^  =  0. 
If  y  is  not  a  limit  point  with  respect  to  the  jth  variable  then  the  solution 
z  e  Rn4^  of  the  linear  system 


a 


(5.14) 


is  a  multiple  of  Ty.  This  provides  a  simple  means  for  the  computation  of  Ty. 

Essentially  all  numerical  continuation  methods  now  are  of  the  predictor- 
corrector  type.  In  other  words,  information  on  an  already  computed  portion 
of  the  curve  is  used  to  calculate  a  suitable  extrapolation  approximating 
a  further  curve-segment.  Some  point  on  this  extrapolating  curve  is  then 
chosen  as  the  starting  point  for  a  corrector-iteration  designed  to  converge 
to  a  point  on  the  solution  curve. 

A  general  system  has  been  developed  for  the  computation  of  the  -solution 
of  (5.8)  through  a  given  starting  point.  It  uses  an  Euler  predictor  and 
Newton's  method  as  corrector.  The  general  scheme  of  the  process  may  be 
described  as  follows: 

1.  [Input]  'starting  point'  y  e  R(F); 

j:=‘ index  such  that  (Ty).  7*0' 


i. 


2.  'From  (5.15)  calculate  tangent  vector'  Ty 

3.  j:=  'index  such  that  j(Ty).|  =  max  | ( Ty ) .  j  * 

J  i*1 . n+1  1 

4.  'predict  steplength'  h 

5.  'compute  Euler  point'  z:=y  +  hTy;  'right  side'  fj; 

6.  'apply  corrector  iteration  to  Fz*c,  (e^)Tz  =  S'; 
normal  return:  y:='final  iterate';  go  to  7; 
error  return:  h:=h/2;  go  to  5; 

7.  vf  'further  step  required'  then  go  to  2 

Let  y^yts.),  i  =  0....,  k,  k  >_  1 ,  be  the  approximations  of  the  desired 
solution  y  =  y(s)  computed  by  this  process.  For  the  steplength  estimation  at 

i/ 

y  we  consider  the  quadratic  interpolant 


(5.15) 


P  ( s )  =  yk  +  (s-sk)Tyk  +  (s-sk)‘'wk 
k  1  /  k-1  k  Tl kv  .  _  ,  |  h  h-1 

w  =  -y  (y  -  y  +  hkiy  )» i»k  -  I  ly  -  y 

hk 


with  p(sk)  =  yk,  p' (sh)  *  Tyk 


)  =  y 


Then  it  can  be  shown  that 


(5.16) 


l y ( s.  +  h)  -  p(s,  +  h)|  L  i  h‘ 


+  0(h  j ,  as  h  -*•  0. 


For  a  given  tolerance  e  >  0 
formula. 

(5.17)  h  =  min 


this  suggest  the  use  of  the  asymptotic  steplength 


In  our  system,  e  is  adaptively  computed  as  c  =  3p  using  a  "safe"  convergence 
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radius  p  of  the  corrector  iteration  and  a  suitable  factor  0  <  0  <  1, 

(see  eg.  [12]).  Figure  5.1  shows  this  schematically.  For  further  details, 
we  refer  to  [11]. 


5.3  A  Numerical  Example 


We  consider  the  boundary  value  problem  (5.1)  with  the  coefficient 
functions 


a(s)  =  yli-  ,  -1  <  s  <  °°  ;  b(s)  =  0,  s  c  R\ 


(5.18) 


f(s)  = 


1  3 

X  for  jj-<s<^-,X>_0, 
0  otherwise 


This  may  be  viewed  as  a  model  of  a  one-dimensional  rod  of  lenght  one  which 
is  clamped  at  the  endpoints  and  subjected  to  a  load  in  the  direction  of  its 
axis.  Then  x  represents  a  measure  of  this  load. 

In  this  special  case,  the  exact  solution  of  (5.1)  for  x  >  0  turns  out  to 


♦  0  1  x  1  4  i 


(5.19)  u ( x )  = 


+  ^  In  (1  +  Xa(x  -  |-)) 


»  4 ix 14 ' 


a(a-ix)x  +  la(l  - - 1 - )  +  1  tl  (1  +  i-a  X  ) ,  r<  X  <  1, 

c  4  1  +  jZ\  A  c  4 


where  a  =  1/(1  -  a)  and  a  =  a(u'(0))  is  the  solution  of 


(5.20) 


(1  +  - Vt-  )  +  ^ln(1  +  =  1* 

1  +  2"  otX 
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Thus  we  have  here 

(5.21)  0^u(x)<^u*(x),  0  <  x  <  1, 

where  the  limiting  function  u*  is  given  by 


(5.22) 


(  ,  0  i  x  <  ~  , 

l  1-x  »  j  <_*  1  1  . 


Now  let  A  be  any  partition  (2.6)  of  I  and  n  =  m(A)-l.  As  in  the  previous 
section  we  introduce  x  as  the  (n+l)st  variable  and  hence  write  (5.2)  in  the 
form  of  the  undetermined  system  (5.8)  (with  c  =  0).  Now  starting  from 
s  =  0,  y=0  the  continuation  process  described  above  was  run  up  to  values  of 
the  load  A  near  30  where  max{u(x),  x  z  1}  >  .8  max{u*(x),  x  e  I}.  Some 
of  the  results  for  uniform  partitions  with  stepsize  h  =  1/12,  h  =  1/24,  and 
h  =  1/48  are  given  in  Tables  5.1  and  5.2  for  the  M-j  2  and  Ml  norm, 
respectively.  The  table  columns  are  as  follows: 

(1)  Parameter  s 

(2)  Load  X 

(3)  Exact  error  | e { ^ 

(4)  Error  estimator  e(A) 

(5)  Relative  error  in  perant,  100  •  |e|.  / | un | , 

I 9p  u  i 

(6)  Efficiency  index  0 


h  =  1/2 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

0 

0 

0 

0 

0 

0 

.5000 

.4880 

. 8336 ( -2 ) 

.8334(-2) 

8.41% 

1 .00 

1.996 

1.955 

. 3502 ( -1  ) 

. 3462 ( -1 ) 

9.39% 

.989 

3.997 

3.936 

. 7339(-l ) 

. 6854 ( - 1 ) 

11.3% 

.934 

6.000 

5.931 

.1077 

.91 02 ( -1 ) 

12.9% 

.846 

12.02 

11.94 

.1695 

.1030 

14.9% 

.608 

18.02 

17.94 

.1936 

.  91 33  ( -1 ) 

15.1% 

.472 

24.02 

23.94 

.2018 

. 7880 ( -1  ) 

14.8% 

.390 

30.02 

29.94 

.2032 

. 6848 ( -1 ) 

14.3% 

.337 

(1) 

(2) 

0 

0 

.4997 

.4767 

1.993 

1.912 

3.994 

3.873 

6.000 

5.863 

12.03 

11.88 

18.04 

17.89 

24.04 

23.89 

30.04 

29.89 

.4072 (-2) 

. 4073 ( -2 ) 

.  1 720 ( -1  ) 

.  1 71 5  ( -1 ) 

.  371 5 ( -1 ) 

. 3645 ( - 1 ) 

.  5695 ( -1  ) 

. 5387 ( -1  ) 

.1014 

.8049( -1 ) 

.1260 

.8359(-l) 

.1395 

.7931 (-1) 

.1469 

. 7332 ( -1  ) 

(5) 

(6) 

0 

0 

4.20% 

1 .00 

4.70% 

.997 

5.78% 

.981 

6.84% 

.946 

8.92% 

.794 

9.86% 

.663 

10.2% 

.568 

10.4% 

.499 

0) 

(2) 

(3) 

(4) 

(5) 

(6) 

0 

0 

0 

0 

0 

0 

.4996 

.4564 

. 1 934 ( -2 ) 

.1949 (-2) 

2.08% 

1.01 

1.987 

1.835 

.8231 (-2) 

.8225( -2) 

2.33% 

1 

.999 

3.986 

3.753 

. 1 809 ( -1 ) 

.1801 (-1 ) 

2.88% 

.995 

5.997 

5.730 

.2846 

. 2804 ( -1 ) 

3.46% 

.985 

12.07 

11.77 

,.5504(-l) 

. 5065 ( -1 ) 

4.86% 

.920 

18.08 

17.78 

.7331  (-1  ) 

.  61 33 ( — 1 ) 

5.74% 

.837 

24.08 

23.78 

.8581 (-1) 

. 6503 ( -1 ) 

6.30% 

.758 

30.08 

29.78 

.  9446 ( -1 ) 

. 6520 ( -1 ) 

6.66% 

.690 

Table  5.2:  "Rod"  Problem,  I.L  0  -norm,  h  =  1/12,  h  =  1/24,  h=l/48 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

0 

0 

0 

0 

0 

0 

.4999 

.4880 

. 251 3 ( -1 ) 

.  2465 ( -1  ) 

19.9% 

.981 

1.996 

1.955 

.1663 

.1488 

31 .8% 

.895 

3.997 

3.936 

.4793 

.3571 

48.9% 

.745 

6.000 

5.931 

.8323 

.5047 

63.3% 

.606 

12.02 

11.94 

1.679 

.5974 

89.7% 

.356 

18.02 

17.94 

2.195 

.5343 

103% 

.243 

24.02 

23.94 

2.527 

.4627 

110% 

.183 

30.02 

29.94 

2.755 

.4029 

114% 

.146 

s 

h  =  1/24 

- -  . 
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\ 

(i) 

(2) 

(3) 

(4) 

(5) 

(6) 

1 

0 

0 

0 

0 

0 

0 

.4998 

.4767 

. 1 237 ( -1  ) 

. 1 226 ( - 1 ) 

10.0% 

.991 

1 

1.993 

1 .912 

. 8539 ( -1 ) 

.8088(-l ) 

16.7% 

.947 

I 

3.994 

3.873 

.2677 

.2308 

27.6% 

.862 

6.000 

5.863 

.5024 

.3860 

38.5% 

.768 

j  j 

12.03 

11.88 

1.184 

.6401 

63.4% 

.541 

18.04 

17.89 

1 .685 

.6806 

78.8% 

.404 

'  1 

24.04 

23.89 

2.045 

.6517 

88.9% 

.319 

30.04 

29.89 

2.310 

.6054 

96.0% 

.262 

h  =  1/48 

• 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

, 

o 

0 

0 

0 

0 

0 

•  1 

.4996 

.4564 

.5906(-2) 

. 5879 ( -2 ) 

5.01% 

.996 

1.987 

1.835 

.4119(-1) 

.  401 5 ( - 1 ) 

8.39% 

.975 

1 

3.986 

3.753 

.1367 

.1273 

14.5% 

.931 

1 

5.997 

5.730 

.2732 

.2397 

21 .2% 

.878 

12.07 

11.77 

.7398 

.5308 

39.8% 

.717 

1 

18.08 

17.78 

1.148 

.6795 

53.8% 

.592 

I 

24.08 

23.78 

1 .478 

.7371 

64.3% 

.499 

1 

30.08 

29.78 

1 .745 

.7478 

72.6% 

.429 

Rod"  Proble 

m,  |.|1>oo  nc 

)rm,  h  =  1/12, 

h=  1/24,  h 

=  1/48 

. i _ 

Lj _ _  _ _ J 
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For  small  X  the  error  decreases  again  linearly  with  h.  For  larger  X 
the  errors  are  still  very  large  and  the  asymptotically  expected  convergence 
has  not  yet  set  in.  As  in  the  linear  case,  the  efficiency  index  improves 
with  decreasing  relative  error  and  0  is  close  to  one  already  for  errors  of 
the  order  of  5%.  Not  unexpectedly,  the  errors  in  the  |.|,  norm  are 
much  worse  than  in  the  l-li  2  norm  ^or  reasonable  errors  the  efficiency 
index  is  again  close  to  one. 

Clearly  the  results  in  these  tables  are  practically  useless  for  higher 
values  of  X .  We  would  need  considerably  smaller  step  sizes  h  to  get 
acceptable  relative  errors  for  large  X.  In  order  to  reduce  the  resulting 
very  heavy  computational  cost,  we  need  to  use  once  again  adaptive  mesh 
constructions. 

For  this  we  apply,  as  in  the  linear  case,  the  principal  of  equilibration 
of  the  error-indicators.  Thus,  for  any  fixed  X  all  those  subintervals  I. 

J 

are  subdivided  for  which  the  error-indicator  exceeds  a  certain  cut-off 
value  n"  .  On  the  new  partition  an  approximate  solution  is  obtained  by 
interpolation  and  then  corrected  by  means  of  a  few  Newton-iteration  steps. 

Now  the  continuation  process  is  applied  again  until  a  further  refinement 
becomes  necessary. 

There  are  as  yet  numerous  open  questions  about  this  combination  of 
continuation  and  adaptive  mesh  refinement.  But  as  Table  5.3  and  Figure  5.2 
show,  the  process  is  highly  effective. 


0 

.5000 

1.999 

3.999 


3.999 

5.996 

7.997 


7.997 

9.996 

12.00 

18.00 


18.00 

20.00 

22.00 

28.00 


28.00 

30.00 


0 

.4959 

1.985 

3.977 


3.977 

5.967 

7.966 


7.966 

9.962 

11.96 

17.96 


17.96 

19.96 

21.96 

27.96 


27.96 

29.96 


1 el 1 ,2 

e(a) 

0 

0 

. 2535 ( -1  ) 

. 2524 ( -1  ) 

.1016 

.9331 (-1) 

.1887 

.1358 

.7371 (-1) 

. 6952 ( - 1  ) 

. 9464 ( -1  ) 

. 8463 ( -1 ) 

.1138 

. 9400 ( -1  ) 

.61 37 ( -1 ) 

. 5658 ( -1 ) 

. 61 88 ( -1 ) 

. 5690 ( -1  ) 

.6297(-l ) 

. 5764 ( -1 ) 

. 6863 ( -1  ) 

. 6052 ( -1 ) 

. 5205 ( -1  ) 

. 4748 ( -1  ) 

.51 77 ( -1 ) 

.471 7(-l ) 

. 51 72 ( — 1 ) 

. 4702 ( -1 ) 

. 5256 ( -1  ) 

.471 1 (-1 ) 

.4244(-l) 

. 3855 ( -1  ) 

.41 84 ( -1 ) 

. 3803 ( -1  ) 

0 

25.2% 

26.9% 

28.8% 


11.3% 

11.3% 

11.7% 


6.32% 

5.81% 

5.53% 

5.36% 


4.07% 

3.95% 

3.89% 

3.75% 


3.02% 

2.95% 


Table  5.3:  "Rod"  Problem,  Adaptive  Mesh  Construction 
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6.  Conclusions 

Although  much  work  still  needs  to  be  done,  we  believe  that  the  approaches 
presented  here  can  provide  satisfactory  answers  to  the  two  questions  (1.1) 
raised  in  the  introduction.  We  summarize  briefly  some  of  the  conclusions 
that  can  be  drawn  from  our  results: 

a)  The  a-posteriori  error  estimates  provide  reliable  information 
for  a  wide  class  of  problems  and  under  a  variety  of  different 
norms. 

b)  In  all  our  experience,  the  estimates  are  acceptable  even  for 
relative  errors  o£  about  5-10%  and  they  improve  with  decreasing 
errors. 

c)  In  the  one-dimensional  case  the  theory  of  the  a-posteriori 
estimates  is  rather  well  established  but  in  the  two-dimensional 
and  even  more  so  in  the  nonlinear  case  there  are  still  a 
number  of  open  theoretical  questions. 

d)  The  error  estimates  provide  a  very  effective  means  for  controlling 
the  adaptive  construction  of  near-optimal  meshes. 

e)  Experience  has  shown  that  the  errors  for  the  adaptively  constructed 
meshes  show  near-optimal  character  even  in  the  presence  of  severe 
singularities. 

f)  Adaptive  mesh  construction  generally  provides  an  overall  savings  in 
computational  cost  and  improves  the  efficiency  index  of  the  error 
estimates. 

g)  For  an  effective  implementation  of  the  adaptive  mesh  construction 
special  attention  needs  to  be  paid  to  all  aspects  of  data 
management.  A  tree  structure  for  the  representation  of  a 
hierarchy  of  meshes  has  proved  itself  to  be  very  advantageous  in 
this  connection  and  allows  also  for  natural  pivoting  by  nested 
dissection. 

h)  Our  experience  with  the  FEARS  project  confirms  that  general-purpose 
finite  element  solvers  can  be  constructed  which  provide  near- 
optimal  solutions  within  a  given  cost  range. 


i 


66 


i)  In  the  design  of  such  finite  element  solvers  it  is  possible  to 
utilize  natural  parallelity  features  of  the  finite  element  method 
and  to  achieve  a  natural  data  segmentation. 

j)  The  a-posteriori  error  estimates  are  also  very  effective  for 
nonlinear  problems  although  the  corresponding  theory  is  as 
yet  far  from  complete. 

k)  For  the  solution  of  nonlinear  equilibrium  problems  of  the  type 
arising,  for  instance,  in  structural  analysis  stable  continuation 
processes  can  be  designed  which  are  not  affected  by  limit  points. 

l)  In  the  case  of  nonlinear  problems,  the  a-posteriori  error  estimation 
and  adaptive  mesh-construction  techniques  can  be  combined  with  the 
continuation  process  with  highly  advantageous  effects  on  the 
computational  cost  and  the  efficiency  indexes  of  the  error  estimates. 
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